# -*- coding: utf-8 -*-
"""
Created on Tue Sep 11 16:52:14 2012

@author: xudi
"""

import numpy as np
import pylab as pl


'''index'''
size=20000
maxtime=60000
M=np.zeros(size)
N=7 #奇数
xl1=1
xl2=3
e0=8.85e-12
u0=1.26e-7
e1=np.zeros(1437)
e2=np.zeros(2489)
e1[:1437]=xl1*e0
e2[:2489]=xl2*e0

'''fdtd index'''

C=3e8
lambc=7.8e-7
dz=lambc/1000
dt=dz/(2*C)
wc=2.417e15
w0=100*wc
H=np.zeros(size-1)
E=np.zeros(size)
J=np.zeros(size)
P=np.zeros(size)
R=np.zeros(maxtime)
T=np.zeros(maxtime)
Tao=wc

'''medium inedx'''
e=np.append(e1,e2)
mm=0
M[:2000]=e0
mm+=2000
for mm in range(2000,size):
    if (mm-2000)%3926==0:
        M[mm:mm+3926]=e
        mm+=3926
        if mm>2000+3926*(N-1)/2-1:
            M[mm:size]=e0
            break
print mm
print M
'''source'''

tw=5e-16# only one can change
t0=3*tw
t = np.arange(0.,1e-13,dt)
ft=np.exp(-(t-t0)**2/tw**2)
x=ft[0:maxtime+1]
y=np.arange(0,maxtime*dt,dt)

'''fdtd'''
qtime=0
for qtime in range(0,maxtime):
    #H[size-1]=H[size-2]    
    H[:size-1]=H[:size-1]-(dt/(u0*dz))*(E[1:size]-E[:size-1])
    H[999]-=x[qtime]*((e0/u0)**0.5)
    J=((1-(Tao*dt/2))/(1+(Tao*dt/2)))*J+(dt/(1+(Tao*dt/2)))*(M*(wc**2)*E-(w0**2)*P)
    E[0]=E[1]
    E[1:size-1]=E[1:size-1]-(dt/(M[1:size-1]*dz))*(H[1:size-1]-H[0:size-2])
    E[1000]+=x[qtime+1]
    
   
    P=P+dt*J
    R[qtime]=E[100]
    T[qtime]=E[19000]
    qtime+=1

'''FFT'''
#sampling_rate=1/dt
#xf1=np.fft.fft(R,maxtime)
#xf2=np.abs(np.fft.fftshift(xf1))
#freqs=np.linspace(0,sampling_rate,maxtime)-sampling_rate/2


#pl.plot(freqs,xf2,'b')
pl.plot(y,T,'r--')
pl.show()

